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Twisting DNA under a constant applied force reveals a thermally activated transition into a state 
with a supercoiled structure known as a plectoneme. Using transition state theory, we predict the 
rate of this plectoneme nucleation to be of order 10 4 Hz. We reconcile this with experiments that 
have measured hopping rates of order 10 Hz by noting that the viscous drag on the bead used to 
manipulate the DNA limits the measured rate. We find that the intrinsic bending caused by disorder 
in the base-pair sequence is important for understanding the free energy barrier that governs the 
transition. Both analytic and numerical methods are used in the calculations. We provide extensive 
details on the numerical methods for simulating the elastic rod model with and without disorder. 



I. INTRODUCTION 

When overtwisted, DNA forms supercoiled structures 
known as plectonemes (as seen in the lower right of 
Fig. [TJ, familiar from phone cords and water hoses. Sin- 
gle molecule experiments commonly hold a molecule of 
DNA under constant tension and twist one end; the ap- 
pearance of a growing plectoneme can be thought of as 
the nucleation of a new phase that can store some of the 
added twists as writhe. Recent experiments that hold 
DNA near this supercoiling transition have shown that it 
is not initiated by a linear instability (as it is in macro- 
scopic objects), but is rather an equilibrium transition 
between two metastable states, with and without a plec- 
toneme. These states are separated by a free energy bar- 
rier that is low enough to allow thermal fluctuations to 
populate the two states, but high enough that the charac- 
teristic rate of hopping is only about 10 Hz. Two exper- 
imental groups, using different methods to manipulate 
the DNA (one, an optical trap [1 ; the other, magnetic 
tweezers [2]), have observed this nucleation at the tran- 
sition and reported similar qualitative and quantitative 
results. 

Understanding the rate of plectoneme nucleation at the 
supercoiling transition is a useful goal for both biology 
and physics. First, the biological function of DNA is 
tied to its microscopic physical characteristics, and plec- 
toneme nucleation is sensitive to many of these. The 
microscopic dynamics of DNA in water is one such fac- 
tor: while often theorized as a cylindrical rod in a viscous 
liquid, these dynamics have not been well studied exper- 
imentally. The nucleation rate is also sensitive to the 
intrinsic bend disorder present in a given DNA sequence, 
potentially providing information to clarify the degree 
of bend disorder, which is debated in current literature 
[3HS1- Popular elastic rod models for DNA can also be 
tested. Second, DNA supercoiling provides to physics a 
unique testing ground for theories of thermal nucleation. 
The theory of thermal nucleation in spatially extended 
systems (critical droplet theory) was essentially proposed 
in its current form by Langer in the 1960's |8 (see 
Hanggi [9[ section IV. F], and also Coleman [10] for the 
corresponding 'instanton' quantum tunneling analogue.) 



Experimental validation of these theories has been diffi- 
cult in bulk systems, however, for reasons that DNA nu- 
cleation neatly bypasses, (a) The nucleation rate in most 
systems is partly determined by the atomic-scale surface 
tension; in DNA the continuum theory describes the en- 
tire nucleation process, (b) Nucleation in bulk phases 
is rare (one event per macroscopic region per quench), 
and hence typically has a high energy barrier. Small es- 
timation errors for this barrier height typically hinder 
quantitative verification of the (theoretically interesting) 
prefactor. In single-molecule DNA experiments, the plec- 
toneme nucleation barrier is only a few and indeed 
rather short segments of DNA exhibit multiple hops over 
the barrier — allowing direct measurements of the tran- 
sition rate in equilibrium, (c) Nucleation in bulk phases 
is normally dominated by disorder (raindrops nucleate 
on dust and salt particles); in DNA the likely dominant 
source of disorder (sequence-dependence) is under the ex- 
perimentalist's control. 

As an illustrative example, consider the classic early 
study of supercurrent decay in thin wires [TTJ [12] . Here 
(a) the superconductor (like DNA) is well described by 
a continuum theory (Ginzburg-Landau theory) because 
the coherence length is large compared to the atomic 
scale, but (c) the rate for a real, inhomogeneous wire 
will strongly depend on, for example, local width fluctu- 
ations. Finally, (b) the rate of nucleation is so strongly 
dependent on experimental conditions that an early cal- 
culation [11 had an error in the prefactor of a factor 
of 10 10 [12] , but nonetheless still provided an acceptable 
agreement with experiment. 

Reaction-rate theory predicts a rate of plectoneme nu- 
cleation related to the energy barrier between the two 
states. We perform a full calculation of this energy bar- 
rier and of the rate prefactor, including hydrodynamics, 
entropic factors, and sequence-dependent intrinsic bend 
disorder, to determine which effects contribute to this 
slow rate. We calculate a rate of order 10 4 Hz, about 
1000 times faster than measured experimentally. The dis- 
crepancy can be attributed to a slow timescale governing 
the dynamics of the measurement apparatus. The exper- 
iments measure the extension by monitoring the position 
of a large bead connected to one end of the DNA, and its 
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FIG. 1: (Color online) The free energy double- well: 
A schematic one-dimensional representation of the two 
met ast able states and the saddle point that separates them. 



II. NUCLEATION RATE CALCULATION 

A. Saddle point energetics 

To understand the dynamics of nucleation, we must 
first find the saddle point DNA configuration that serves 
as the barrier between the stretched state and the plec- 
tonemic state. We model the DNA as an inextensible 
elastic rod, with total elastic energy 

Elastic = ds[^(s) 2 + jT( S )% (1) 



dynamics are much slower than that of the DNA strand 
— thus the bead hopping rate that the experiments mea- 
sure is much slower than the plectoneme nucleation hop- 
ping rate that we calculate. Future experiments may be 
able to slow the plectoneme nucleation rate enough that 
the bead dynamics is unimportant, such that the bead 
motion would reveal the dynamic characteristics of the 
DNA itself. 



We begin in section [II] with the nucleation rate calcu- 



lation in the absence of disorder: |II A| gives the saddle- 
point energy, II B gives the technique for calculating the 
prefactor, |II C| overviews the dynamics of DNA in water, 
and |IID| gives the transition-state theory calculation in 
full. Section III presents the results of the undisordered 
calculation, with qualitative explanations of the magni- 
tudes of the various terms, and shows that the results are 
incompatible with the experiments. This motivates our 
discussion of base-pair disorder in section |IV[ where we 
estimate the disorder-renormalization of the elastic con- 
stants in sectio n |IV B , and formulate and calculate the 
rates in sections II V Cl and flV D[ We conclude in sectionlVl 



where s is arclength along the rod, (3 and T are the lo- 
cal bend and twist deformation angles, respectively, and 
L is the contour length of the rod (for more details, see 
Appendix [Aj the d ynam ical equations of motion will be 
discussed in section II C ) . Fain and coworkers used varia- 



tional techniques to characterize the extrema of this elas- 
tic energy functional [27 ; this revealed a "soliton-like 
excitation" as the lowest-energy solution with nonzero 
writhe [13]. They found that the soliton's energy dif- 
fered from that of the straight state by a finite amount 
in the infinite- length limit. This soliton state, depicted 
at the top of the barrier in Fig. [I] is the one we identify 
as the saddle configuration. 

The shape of the saddle configuration is controlled by 
the bend and twist elastic constants B and C, as well as 
the torque r and force F applied as boundary conditions. 
Defining the lengths 
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we can write the Euler angles characterizing the saddle 
configuration as [T3] [28] 



We draw the reader's attention in particular to the ap- 
pendices, where substantive, general-purpose results are 
presented for numerical discretization and calculations 
with the elastic rod model. Appendix [B] reformulates 
the Euler-angle description in terms of more geometri- 
cally natural rotation matrices. Appendix [C] explains 
how to transform a DNA with N segments from the 37V- 
dimensional Euler-angle or rotation-matrix space to the 
47V-dimensional (x, y, z, <p) space of the natural dynam- 
ics, and the Jacobians needed to transform path integrals 
over the latter into path integrals over the former. Ap- 
pendix [D] discusses the discretization and the rotation- 
invariant forms for bend and twist in terms of rotation 
matrices. Finally, Appendix [E] provides our numerical 
implementation of randomly-bent DNA, mimicking the 
effects of a random base-pair sequence. 
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where s is arclength along the DNA backbone. 

For the experiments in which we are interested, 
L/t ^> 1 (the soliton "bump" is much smaller than 
the length of the DNA), such that we can safely remove 
the soliton from the infinite length solution and still have 
the correct boundary conditions [0(±L/2) = 0, such that 
the tangent vector t points along the z axis at the ends of 
the DNA]. In this case, the linking number in the saddle 
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state is given by [T3] 
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where we have used Eq. ([5| and separated the linking 
number into twist (the first term) and writhe: 
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To find the saddle configuration's torque at the super- 
coiling transition [29 j, we numerically solve Eq. (|6| for r 
using the experimentally observed critical linking number 

k* s m- 

The energy barrier AE is the difference in elastic en- 
ergy between the saddle and straight states at the same 
linking number K s [31]. We find [32] 
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Inserting the experimental values listed in Table [I] into 
Eq. (l8|), we calculate an energy barrier 



AE = 5.5 kT. 



(9) 



This barrier would seem surprisingly small considering 
that typical atomic rates are on the order of 10 13 Hz: 
using this for the attempt frequency in an activated rate 
would give 10 13 e -5,5 = 10 11 Hz for the hopping rate. The 
next sections present a more careful calculation, which 
shows that the timescale for motion over the barrier is 
in fact many orders of magnitude smaller than 10 13 Hz 
due to the larger length scales involved (and even smaller 
when we calculate the bead hopping rate), but also that 
the entropy from multiple available nucleation sites sig- 
nificantly lowers the barrier. 



B. Transition state theory: the basic idea 

When, as in our case, the energy barrier is much larger 
than the thermal energy the rate of nucleation is sup- 
pressed by the Arrhenius factor exp (— AE/kT). Going 
beyond this temperature dependence to an estimate of 
the full rate, however, requires a more detailed calcu- 
lation. We will follow the prescription from Kramers' 
spatial-diffusion-limited reaction-rate theory [9] to cal- 
culate the rate of hopping. The requirements are that 
(1) the timescales involved in motion within the two 
metastable wells are much faster than the timescale 
of hopping, and (2) (for Kramers' "spatial-diffusion- 
limited" theory) the system is overdamped, in the sense 
that the ratio of damping strength to the rate of un- 
damped motion over the barrier top is large. We check 
that these requirements are met for the intrinsic DNA 
nucleation rate after the calculation in section III D[ 



Under these two conditions, Kramers' reaction-rate 
theory tells us that the rate of hopping over the barrier is 
controlled only by the rate of motion through the "nar- 
row pass" at the top of the barrier, since it is much slower 
than any other timescale in the system. This means that 
the hopping rate should be the characteristic rate of mo- 
tion across the barrier top times the probability of finding 
the system near the barrier top, which, in terms of the 
curvature in the unstable direction away from the saddle 
point, we can write schematically as 

^hopping = (Characteristic rate of motion at barrier top) 



x (Prob. of being at top) 



(10) 



/Energy curvature \ . 

= — (Prob. ol being at top). 

\ Damping J 

(ii) 

It is important to note that, in current experiments, 
the measuring apparatus violates condition (1) above. 
The measurement of the extension is only an indirect 
readout of the configurational state of the DNA — it is 
a measure of the position of a large bead connected to 
one end of the DNA strand. If the bead has much slower 
dynamics than the DNA, then it will set the ch aracter- 
istic rate of motion in Eq. (10). In section IIID 1 we will 



find that this is the case for the experimental numbers 
we use. Therefore, the rate we will calculate is a plec- 
toneme nucleation hopping rate that is not the same as 
the (slower) bead hopping rate. We will find also that fu- 
ture experiments may be able to measure the underlying 
plectoneme nucleation hopping rate by testing regimes 
where the bead hopping rate is not limited by the bead 
dynamics. 



C. Dynamics of DNA in water: the diffusion tensor 

To find the rate of motion over the barrier top, we need 
to know the microscopic dynamics. We will be treat- 
ing the DNA strand as a series of cylindrical segments, 
parametrized by the Cartesian coordinates (x n: y n ,z n ) 
of one end of each segment plus the Euler angle ip n 
that controls local twist (see section CI for a discus- 
sion about the choice of coordinates). Assuming over- 
damped motion such that we can neglect inertial terms, 
we will write the equations of motion in the form [with 



r n = (x n , y n , z n , il> n ), i,j labeling coordinates, and m, n 
labeling segments] [9] 



dt 



-Mo 



mz,nj 



dE 

dv n j 



(12) 



M is the diffusion tensor, which transforms forces to ve- 
locities. 

The simplest diffusion tensor produces motion propor- 
tional to the local forces, making it diagonal in segment 
number n and coordinate i: 



ti ^-diagonal 
mi,nj 



J_ A 
dX U ' mn 



for i,j e {x,y,z} 
for i = j — ip 



(13) 
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TABLE I: Parameter values for nucleation rate calculation. 



Symbol 


Description 


Value 


TD 
J J 


bend elastic constant 


( a q ^ w~r m 
^4o nmj Kl m 


O 


twist elastic constant 


(oy nmj /ci yj 


TP 

r 


applied force 


i.yo piN yj 


kl 


thermal energy at 23.5 C 


4.09 pN nm 


T 

L 


basepair length of DNA strand 


f 4U nm yj 


If* 


critical linking number 


8.7 [I] 


r 


saddle point torque [Eq. JoftJ 


25 pN nm 


o 
I 


soliton length scale [Eq. H] a 


13 nm 


TD 
1\ 


bead radius 


zou nm 


T] 


viscosity of water at 23.5° C 


Q 99 V 10 — 10 t^1\T a / nm 2 

y.zz x iu piN s / nm 


N 


number of segments 


740 


d 


length of segment 


1 nm 




DNA hydrodynamic radius 






1.2 nm 


c 


translation viscosity coeff. [Eq. ( 


14 


)] 


1.54 x 10" 9 pN s / nm 2 


A 


rotation viscosity coeff. [Eq. (|15 


1] 


1.67 x 10~ 8 pN s 



a These values have been calculated using B m = B, that is, for disorder D — 0. See Eq. (30) 



The viscous diffusion constants are set so they reproduce 
the known diffusion constant for a straight cylinder of 
length B/kT (the bending persistence length) and radius 
r D = 1.2 nm: [THUS] 



s \n(B/(kTr D )) v 7 

A = 2ir7ir 2 D , (15) 

where the viscosity of water rj — 9.22 x 10 -10 pNs/nm 2 
at the experimental temperature of 23.5° C. 

It is important to note that the DNA in this exper- 
iment is attached to a large bead (with radius R « 
250 nm) that must also be pulled through the water dur- 



ing the transition [33]. We take this into account by 
setting the translational diffusion constant for the final 
segment in the chain according to Stokes' Law: 

M NijNj = 5ij/(67rr}R) for i,j G {x,y,z}. (16) 

Hydrodynamic effects may also be important, which 
introduce interactions between segments: as segments 
move, they change the velocity of the water around them, 
and this change propagates to change the viscous force 
felt by other nearby segments. Following Ref. [16], we in- 
corporate hydrodynamic effects by using a Rotne-Prager 
tensor for the translational diffusion, modeling the strand 
as a string of beads [34] : 
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32 a 



> 2a, m 7^ n 



3 r n 



Dq Si 



l 3 ^ 32 a 

for r mn < 2a, m ^ n 
for m = n 



3^ 



(17) 



for i : j G {x : y,z}, where Do = (67r^a) _1 and a is an effective bead radius chosen such that a straight configuration 
of Kuhn length Lk = 2 B/kT (with a number of beads L/Lk) has the same total diffusion constant as a cylinder of 
length L and radius r& (see [H]). With the parameters in Table |IJ we use a = 0.98 nm. 



D. Transition state theory: full calculation 



In the full multidimensional space inhabited by our model, the saddle configuration will have a single unstable 
direction that locally defines the "reaction coordinate" depicted in Fig. [I] The direction of the unstable mode can be 
found numerically by locally solving the equations of motion Eq. ( [12] ). First, the local quadratic approximation to 
the energy is provided by the Hessian 
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d 2 E 



(18) 



where the derivatives are taken with respect to the unit- 
less variables r = x/£o,y/£o : z/£q : i/j (wh ere l p is an ar- 
bitrary length scale [35]; see also section [UT] ). Inserting 
the quadratic form defined by H at the saddle point into 
Eq. (12), we then diagonalize the matrix MH sa dd\e to 



find the dynamical normal modes of the system; the sin- 
gle mode u with a negative eigenvalue — A& is the unstable 
mode at the top of the barrier: 



(19) 



and A5 defines the characteristic rate of Eq. (11). We 



have checked that we find the correct saddle configura- 
tion and unstable mode u by perturbing forward and 
backward along u and numerically integrating the dy- 
namics of Eq. ( 12 ) — one case ends in the straight state 
well and the other in the plectonemic state well. This 
generates the transition path connecting the two wells, 
as illustrated in Fig. |2j 

To find the probability of being near the top of the bar- 
rier in the multidimensional case, we need to know not 
only the energy barrier AE, but also the entropic factors 



coming from the amount of narrowing in directions trans- 
verse to the transition path [36 , which are controlled 
by the remaining eigenvalues of i^ sa ddie an d the Hessian 
^straight of the straight state. The full result from spatial- 
diffusion limited multidimensional transition state theory 
is (to lowest order in kT) [9] 



7 _ Afr / detilf straight /(27rfcr) _ AE / } 

topping - 2?r y | detffsaddle/(27rA;T) | e 

or, in terms of the eigenvalues of each Hessian, 



"^hopping 



2tt' 



n 



4(7V- 
1 



2) ^straight 



/(2tt/cT) 



Mn 



4(7V- 
1 



-AE/kT 



1 Af addle /(27rfeT)| 



(21) 

These are correct if each eigenvalue is sufficiently large 
such that the local quadratic form is a good approxima- 
tion where E < kT. 

In our case, we must deal separately with the two zero 
modes due to invariance with respect to location s s and 
rotation angle p s of the saddle configuration's bump. Ex- 
tracting these directions from the saddle integral, we have 



hopping 



2tt 



T T dSs J 



n 



1 



2tt \ 27T £ JsJp J 2-KkT 



4(iV- 2 ) A st raig ht /(27rfcr) 



\ int3" 2) Ar ddle /(2^r)i 



det -^straight ^ — AE/kT 
I det' Saddle | 



(22) 



r 



where the Jacobians J s = \df s /ds s \ and J p = \dr 3 /dp 3 \, 
and det 7 represents the determinant without the two 
zero modes (but including the unstable mode). Numeri- 
cally, J s and J p are calculated using the known forms for 
derivatives of the saddle point's Euler angles a s [Eqs. ([5|] 
with respect to s s and p s : J s = \ \J T (f. g )] ~ 1 da s / ds s \ an d 
J p — \ [J T (f s )]~ 1 da s /dp s \, where J is defined in Eq. ( |C5[ ). 

We can now check that we meet the requirements for 
using Kramers' theory set out in section II B First, we 



check condition (1) by looking at the smallest nonneg- 
ative eigenvalues of MH. The slowest mode is trans- 
verse motion with wavelength 2L, and since the bead 
has much larger viscous drag than the rest of the DNA 
chain, it sets the damping for this motion; this produces 
a frequency F/(67rr]RL) w 600 Hz. The other modes all 
have frequencies of order 10 4 Hz or faster. We will find 
that the calculated Kramers rate of hopping lies between 
these two timescales — this means that, while the bead 
motion is too slow to follow the fast hopping, Kramers' 
theory should correctly give the plectoneme nucleation 
hopping rate for a fixed bead position. We can check 



(2) by comparing the characteristic rate for undamped 
barrier motion to the characteristic damping rate. In the 
spirit of section IIIB[ we are dealing with a portion of 
DNA of length £3, such that it has a mass pis, where 
p = 3.3 x 10 -21 g/nm is the linear mass density of DNA 
[T7] , and the energy curvature at the barrier top is on 
the order of tt 4 B/£^. Then the damping coefficient is 
(/p = 5 x 10 11 s _1 , and the rate for undamped bar- 
rier motion is ^(7r 4 B/£ 3 B )/(p£ B ) = 4 x 10 8 s _1 . Since 
the ratio of these values is much greater than one, we 
are firmly in the over damped regime, and Kramers' rate 
theory applies [9]. 
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FIG. 2: (Color online) Snapshots along the transition path. 
After perturbing the saddle state along the unstable direc- 
tion u (bottom left and top right snapshots), we integrate 
the equations of motion in Eq. ( 12 ) to follow the unstable 
dynamics into the two metastable wells. Here, we use Rotne- 
Prager dynamics with D — 0, the timestep between frames is 
1.5 x 10 -6 s, and other parameters are as given in Table [i] 



FIG. 3: (Color online) Transition state theory plectoneme 
nucleation hopping rate (bottom) and the factors that con- 
tribute to that rate (top) versus external force. (Blue; lower) 
Energetic factor exp (— Ai^/ZcT); (Red; middle) entropic fac- 
tor exp(Sy/c); (Green; upper) dynamic factor Afe/27r, in Hz; 
(Black) the final hopping rate per unit length from Eq. (22). 
The calculation is performed for the experimental conditions 
in Ref. [T], with L = 2.2 kbp (squares) and L = 4.2 kbp (di- 
amonds). On this log scale, adding the three distances from 
the horizontal line at 10° in the top plot produces the final 
rate. Note that the entropic factor cancels the slowing from 
the energy barrier factor. (Here the calculation is performed 
without intrinsic bend disorder, produci ng s mall or even neg- 
ative free energy barriers — see section 



IV 



III. 



INITIAL RESULTS AND ORDER OF 
MAGNITUDE CHECKS 

A. Initial results 



We calculate the rate in Eq. (22 ) using numerical meth- 



ods described in the Appendices. We will quote the re- 
sults of the calculation by looking individually at the fac- 
tors that contribute to the rate. Writing the rate in var- 
ious simple forms, 



u . i± -AE/kT+S/k 

^hopping — c 
Z7T 

^ -AF/kT. 



2tt' 



(23) 
(24) 



S encapsulates the entropic factors coming from fluctu- 
ations in the straight and saddle configurations, and the 
effective free energy barrier AJ 7 = AE — TS provides the 
relative probability of being near the top of the saddle [9] . 



For the experimental parameters in Table [IJ using 
Rotne-Prager dynamics [37], we find Xb/^Tr = 4.0 x 
10 4 Hz, AE/kT = 5.5, and S/k = 5.8, such that 
AT jkT = — 0.3. There are two surprises here: (1) 
the characteristic rate of motion over the barrier is very 
slow compared to typical atomic timescales, and (2) the 
entropic factors are so large that they completely erase 
the energy barrier. We consider these issues in more de- 
tail in the next three sections. We will find that the slow 
characteristic rate comes from the larger length scales 
involved in the transition, and that entropy wins over 
energy due to the length of the DNA. (Roughly speak- 
ing, since we calculate a probability per unit length of a 
plectoneme critical nucleus, for long enough DNA there 
will always be such a nucleus.) 

These rate factors and the corresponding rate per unit 
length of DNA are plotted as a function of external force 
for the two experimentally-measured lengths in Fig. [3] 
[38] . This rate per unit length itself depends on the 
length because (1) the saddle-state torque r changes with 
L and (2) the twisting component of the unstable mode 
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FIG. 4: (Color online) Unstable mode at top of barrier, (top) 
The four components of the unstable mode eigenvector as 
a function of arclength s along the DNA strand, (bottom) 
Plot of dr — yj 'dx 2 + dy 2 + dz 2 ] the peaks show the locations 
where the (Cartesian) motion of the DNA is greatest when 
traversing the barrier. Note that the width of the pe ak is 
about 75 nm — inserting this length scale into Eq. ( 25 ) pro- 
duces a prefactor ~ 10 5 Hz. 



difj is length-dependent at the lengths of interest 



B. Order of magnitude estimates of the dynamical 
prefactor 

Typically, rates for atomic scale systems have prefac- 
tors on the order of 10 13 Hz. Indeed, inserting typical 
atomic length scales (A) and energy scales (eV) into the 
simple rate equation Eq. ( pT| ), and using Stokes' law for 
an angstrom sized sphere in water, the energy curvature 
is 1 eV/A 2 , and the damping is 67rr]r a = 10 -13 eV s/A 2 , 
producing a dynamical prefactor of 10 13 Hz. 

Why, then, is our dynamical prefactor of order 10 4 — 
10 5 Hz? It turns out that the relevant length and energy 
scales for the DNA supercoiling transition are not atomic. 
For the saddle state, we are dealing with length scales on 
the order of tens of nanometers (much larger than single 
atoms), and energy scales related to the elastic constants 
(B/£ « 10 pN nm < 0.1 eV). 

To arrive at a better estimate, we can approximate 



the saddle energetics from bending energies only. Con- 
sider approximating the saddle state as a straight con- 
figuration with a single planar sinusoidal bump of length 
£b and amplitude A. Since the elastic bending energy 
is Eb = y ) 2 ^ s > where the relevant component 
of the tangent vector is in this case t — Ak t cos (k t s) 
for wavenumber k t = tt/jHb, the total bending energy 



Bt B ju4 \ 2 
2 K t /1 



tt 4 B A 2 . 



this leads 



for the bump is Eb — 2 — ^ 2 

to an energy curvature with respect to amplitude of 
d 2 E B /dA 2 = 2LgB. The viscous damping coefficient cor- 

B 

responding to a rod of radius tb and length £b moving 
sideways through water is £b(, with £ given by Eq. (14). 



Putting this together, our back-of-the-envelope estimate 
for the prefactor is [see Eq. (11)] 



Energy curvature ir 4 B 



Damping 



4c" 



(25) 



We see that the prefactor is strongly dependent on the 
length scale ^ of the bending of DNA in the unstable 
mode motion. This length scale should be related to the 
length £ ~ 10 nm, defined in Eq. Q, characterizing the 
shape of the saddle configuration. As shown in Fig. [4| 
we can check the amount of DNA involved in the unsta- 
ble mode motion by looking at the unstable mode eigen- 
vector. This reveals that a better estimate for £b is in 
fact 75 nm; inserting this into Eq. ( 25 ) gives a prefactor 
~ 10 5 Hz, agreeing with the order of magnitude found in 
the full calculation. This simple calculation shows, then, 
how the 8 orders of magnitude separating the atomic 
scale rates from that of our full DNA calculation arise 
from the smaller energy scales and larger length scales 
involved. 



C. Understanding the entropic factor 

We calculate an entropy S that entirely cancels the 
energy barrier AE. What sets the size of S7 Comparing 
Eq. (22) and Eq. (23), we see that the entropic factor 
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(26) 



This factor comes from comparing the size of fluctua- 
tions in the normal modes of the straight and saddle 
states. We expect that most of the modes will be sim- 
ilarly constrained in the two states, except for the two 
zero modes that appear in the saddle state. These zero 
modes create a family of equivalent saddle points at dif- 
ferent locations and rotations along the DNA, each of 
which contributes to the final rate. Imagining counting 
the number of equivalent saddle points along 2tt radians 
in p and L nanometers in 8, we can write the entropic 
factor in the form 



o(S/k) 



2tt L 

Po So' 



(27) 
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FIG. 5: (Color online) Bound on free energy barrier from 
experimental extension distribution, (top) A histogram of 
measured extensions near the supercoiling transition, for F — 
2 pN, L — 2.2 kbp [HUH], clearly demonstrating bistability. 
(bottom) Fitting the negative natural logarithm of the proba- 
bility density to three quadratic functions indicates that there 
is a free energy barrier separating the straight state (longer 
extensions) from the supercoiled state (shorter extensions) of 
at least 2 kT. 



where po and so define how far one must move the soliton 
bump along p and s to get to an independent saddle 
point [40 j. We expect that po should be about tt radians 
(giving two independent saddle points at each s), and sq 
should be of the order of the length scale of the soliton, 
£ w 10 nm, producing pqSq ~ 30 nm. And indeed, using 
S/k = 5.8 found in the full calculation, Eq. (27) gives 



PoSo = 14 nm [41]. Thus the size of the entropic factor 
makes sense: it is large because there are many equivalent 
locations along the DNA where the plectoneme can form. 



D. Estimates of the free energy barrier and bead 
dynamics 

The effective free energy barrier in our calculation is 
reduced to near zero. This is due to the entropic fac- 
tor, which favors the saddle state due to the location 
and rotation zero modes [42 . With a free energy bar- 
rier this small, though, plectonemes would form sponta- 
neously even at zero temperature, and with no barrier to 
nucleation, no bistability would be observed — indeed, 
this would violate our original assumptions necessary for 
the use of transition state theory itself. However, the 
fact that bistability is observed in experiment [TJ [2] as- 
sures us that the effective free energy barrier is in reality 



nonzero; furthermore, the degree of bistability can give 
us a reasonable bound on the size of the barrier. 

Two separate experimental groups have directly mea- 
sured the distribution of extensions observed for many 
seconds near the supercoiling transition; one such his- 
togram [TJ [18] is shown in the top of Fig. [5j and the 
distributions measured by the other group [2]appear re- 
markably similar. Taking the natural logarithm of this 
probability density produces an effective free energy land- 
scape in units of /cT, shown in the bottom of Fig. [5] 

To compare the free energy barrier apparent from the 
extension data to the one from our calculation [defined in 
Eq. (23)], there are two subtleties to consider. First, the 
measurements with extension near the middle, between 
the two peaks, are not guaranteed to correspond to con- 
figurations that are traversing the saddle between the 
two wells (that is to say, extension is not the true reac- 
tion coordinate). Since adding extra probability density 
unrelated to the transition near the saddle point would 
lower the measured effective free energy barrier, we will 
only be able to put a lower bound on the true AF. Sec- 
ond, looking at the transition state rate equation for one- 
dimensional dynamics [9], 



feiD 
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Veil -AE/kT 
^saddle 



(28) 



we see the entropic factor coming from the ratio of energy 
curvatures in the saddle and well states [43]. Comparing 



this to Eq. ( 24 ) , we see that the comparable effective free 



energy barrier should be corrected by this entropic ratio, 
such that 



AF/kT = AE/kT - - In _^L, (29) 

2 A sa dcli e 



where the As are the curvatures in the well and saddle 
states. As shown in the bottom of Fig. [5j we can use 
fit parabolas to estimate this entropic correction, finding 
0.5 kT. Using the heights of the parabolic fits, we find 
that AE = 2.5 fcT, such that our lower bound on the 
effective free energy barrier is about 2 kT. 

Finally, we can now check whether the bead dynamics 
slow the hopping measured in experiments. The curva- 
ture of the parabola at the top of the barrier in Fig. [5] 
gives A sa ddie/27r = 3 x 10 -3 pN/nm; this matches with 
the energy curvature F/L that controls the bead's mo- 
Thus, as in section IIID 



tion in section 



ira 



the char- 
acteristic rate of bead motion is about 600 Hz — using 
the lower bound of 2 kT for the free energy barrier then 
gives a bead hopping rate of 80 Hz, near the observed 
hopping rate. This provides an explanation for the dis- 
crepancy between the fast rate we calculate and the slow 
measured rate: the experimental rate is limited by the 
dynamics of the bead. 
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IV. INCLUDING INTRINSIC BENDS 

DNA with a random basepair sequence is not perfectly 
straight, but has intrinsic bends coming from the slightly 
different preferred bond angles for each basepair. This 
can profoundly affect our calculation by both providing 
pinning sites for plectonemes and changing the relevant 
effective viscosity. 



A. History of intrinsic bend measurements 

Since thermal fluctuations also bend DNA, the degree 
of intrinsic bend disorder is difficult to measure, but can 
be estimated using specific DNA sequences that are in- 
trinsically nearly straight. The contribution to the bend 
persistence length from quenched disorder alone (P) can 
be found using the relation [3 (B/kT) -1 — = 
P" 1 + P _1 : P m is found using an intrinsically straight 
sequence (such that P _1 = 0) and compared to P e ff from 
a random sequence. 

Using estimates of wedge angles along with sequence 
information, Trifonov et al. estimated P = 216 nm [4j|5]. 
An experiment using cryo-electron microscopy [3] found 
P m « 80 nm and P e ff « 45 nm, giving an intrinsic bend 
persistence length of P ~ 130 nm. More recently, a 
group using cyclization efficiency measurements found 
P rn — 49.5 ± 1 nm and P e ff = 48 ± 1 nm, from which 
they conclude that P > 1000 nm [6], in striking contrast 
with the previous estimates. 

We include intrinsic bend disorder in our simulations 
by shifting the zero of bending energy for each segment by 
a random amount, para meter izing the disorder strength 
by D = P~ 2 (see section [eT| for a detailed description). 
We are able to locate the new saddle point including dis- 
order, as illustrated in Fig. [6| using numerical methods 
described in section |E3| Due to the disagreement in the 
literature about the correct value of P, we treat it as 
an adjustable parameter and examine the effects of dis- 
order in a range from P = 1000 nm to P = 130 nm 
(D = 0.03 nm" 1 / 2 to D = 0.09 nm" 1 / 2 ). 







FIG. 6: The saddle state with increasing intrinsic bend 
disorder. (Color online) From top to bottom, D = 
(0,0.02,0.04,0.06,0.088) nm" 1/2 , corresponding to persis- 
tence lengths P = (oo, 2500, 625, 278, 130) nm, respectively. 
The saddle state is located numerically by searching locally 
for zero force solutions. 



elastic parameters 
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B. Renormalization of DNA elastic parameters 

It is important to note that the measured elastic con- 
stants B and C are effective parameters that have been 
renormalized by both thermal fluctuations and intrinsic 
bend disorder. Our simulations do not explicitly include 
thermal fluctuations, but incorporate them by using the 
measured effective elastic constants. When we explic- 
itly include bend disorder, however, we must use micro- 
scopic constants B m and C m adjusted so they create the 
same large-scale (measured) effective constants. Nelson 
has characterized the first-order effect of disorder on the 
elastic constants [5] ; correspondingly, we use microscopic 



C. Rate equation with disorder 

If the disorder is large enough, plectoneme formation 
will be strongly pinned to one or more locations along the 
DNA. In this case, the zero modes have vanished, and we 
find the total rate by adding the contributions from each 
saddle point at each location Sj. Using Eq. (20), 



^hoppin£ 



\,j / det Straight e ~(E Si 
27T V | det Saddle I 



t)/kT 



(32) 



We can determine when this approximation will be 
valid by checking that fluctuations in the saddle point 
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position Sj are small compared to the spacing between 
locations. The size of fluctuations in Sj can be found in 
the quadratic approximation as 



ASj = JkT / ^2^saddle(5j) 



j kT / _f_ ( D ^saddle (Sj) 

I ds 2 \ dD 



(33) 



where we have replaced Caddie by its f irst- order approx- 
imation at low disorder D (see section E2). Calculating 



the seco nd d erivative numerically at the pinning sites us- 
ing Eq. |e8|, we find £2 ^^^ ~ 0.5 pN/nm. We 

can thus avoid special treatment of the translation modes 
when the fluctuations in sj are much smaller than the av- 
erage distance between pinning sites, about 75 nm (see 
Fig. 14). Setting Asj < 75 nm in Eq. (33) then produces 
a lower bound on the disorder strength D (or, equiv- 
alently, an upper bound on the intrinsic bend disorder 
persistence length P): 



D > 10" 3 nm" 1/2 , or P < 10 6 nm. 



(34) 



The experimental estimates for the disorder persistence 
length P are typically much smaller than this bound 
(hundreds to thousands of nanometers), so the large dis- 
order limit [Eq. (|32|)] should be valid for our calculation. 



D. Results with disorder 

Fig. [7] displays our results for the hopping rate and ef- 
fective free energy barrier as a function of the intrinsic 
bend disorder strength D, for D in the range correspond- 
ing to experimental estimates of the persistence length P. 
The hopping rate is calculated using Eq. ( [22] ) at zero dis- 
order and Eq. (32) with disorder (in this case summing 



over 10 saddle points — see section |E3| for details). The 
effective free energy barrier is calculated using Eq. ( [24] ) 
(using the average A5 over the 10 saddle locations). We 
see first that the nucleation rate is not significantly al- 
tered by the intrinsic disorder, remaining between 10 4 
and 10 5 Hz (too fast to measure with current experi- 
ments). The effective free energy barrier, however, rises 
above zero with increasing disorder, making our calcula- 
tion more physically plausible. Furthermore, note that 
only the larger experimental estimate for intrinsic disor- 
der (P = 130 nm; green dashed line) is consistent with 
our lower bound on AJ 7 of 2 kT. 

In Fig. |8j we plot the components that contribute to 
the rate, defined in Eq. ( [23] ) : the dynamic factor A& 
(green), the entropic factor exp(S/k) (red), and the en- 
ergetic factor exp (—AE/kT) (blue). To explore the vari- 
ance caused by sequence dependence, we calculate these 
factors for five different random intrinsic bend sequences, 
for a single plectoneme location s (here, the location s* 
predicted to have the lowest energy barrier by first-order 
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FIG. 7: Hopping rate and effective free energy barrier vs. dis- 
order magnitude, for conditions in Table [I] (Color online) The 
green dashed line corresponds to one experimental estimate of 
the bending order persistence length, P — 130 nm [3], and the 
blue hatched region corresponds to another, P > 1000 nm [6]. 
The hopping rate ^hopping (top) does not change significantly 
with the addition of intrinsic bend disorder. The effective free 
energy barrier A J 7 (bottom) increases with disorder, such that 
only the smaller P is consistent with the lower bound of 2 kT 
found using the distribution in Fig. [5] The points show re- 
sults from a single random DNA sequence, and the error bars 
are estimates that include uncertainty in DNA sequence and 
elastic constants, corresponding to the range of values found 
in Fig. |8| 



perturbation theory; see section E2). Depending on the 



actual degree of disorder, sequence dependence creates a 
spread in the hopping rate of around one order of mag- 
nitude. 

Since the hopping rate is exponentially sensitive to en- 
ergy scales at the transition, it will also be important 
to carefully consider our knowledge of the true elastic 
constants B and C. Our values [ B = (43 ± 3 nm)/cT 
and C = (89 ± 3 nm)/cT ] were obtained directly from 
the same experimental setup that produced the hopping 
data PQ, and come from fitting force-extension data to 
the worm-like chain model [19] . The uncertainties in pa- 
rameters correspond to ranges of rate predictions — we 
numerically check these ranges by performing the rate 
calculation using both the upper and lower limits of the 
ranges for the quoted value of the two elastic constants. 
We find that changes in the elastic constants mainly af- 
fect the rate through the energy barrier AE [Eq. Q], 
which is much more sensitive to B than to C. The hori- 
zontal bars in Fig. [5] show the results of changing B from 
its lower to its upper limit; we see that the uncertainty in 
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FIG. 8: (Color online) Hopping rate factors versus disorder 
magnitude for different sequences (for the one best plectoneme 
location s* when D > 0, and with F = 2 pN, L = 2.2 kbp). 
Colors are the same as in Fig. [3] Different solid markers cor- 
respond to five different random seeds (five different basepair 
sequences). Horizontal bars correspond to varying the bend 
elastic constant B for one of the sequences by adding and sub- 
tracting the uncertainty in its measurement, (3 nm) kT. On 
this log scale, adding the three distances from the horizontal 
line at 10° produces the final contribution to the rate from 
location s*. (Including multiple plectoneme locations further 
increases the entropic factor.) Both sequence dependence and 
parameter uncertainty increase the spread of possible rates by 
about an order of magnitude. 



the bending elastic constant produces variations on the 
same scale as the sequence dependence. 



V. DISCUSSION AND CONCLUSIONS 

To calculate the rate for plectoneme nucleation at the 
supercoiling transition, we first use an elastic rod the- 
ory to characterize the saddle state corresponding to the 
barrier to hopping. Using reaction rate theory, we then 
calculate the rate prefactor, including entropic factors 
and hydrodynamic effects. We also analyze the effect of 
intrinsic bend disorder, which simultaneously lowers the 
energy barrier and increases the entropic barrier. We 
find that the experimental rate is in fact set by the slow 
timescale provided by the bead used to manipulate the 
DNA, with an intrinsic plectoneme nucleation hopping 
rate about 1000 times faster than the measured bead 
hopping rate. 

Further insight is gained by studying the factors that 
contribute to the plectoneme nucleation rate. First, the 
energy barrier is calculated analytically using elastic the- 
ory [Eq. Second, the rate of motion at the barrier 
top can be obtained in the full calculation, and the order 
of magnitude (10 5 Hz) agrees with the expected rate of 
motion of a rod in a viscous fluid when inserting the ap- 



propriate length and energy scales. Third, the entropic 
contribution to the prefactor significantly lowers the free 
energy barrier, in a way directly related to the saddle con- 
figuration's translational zero mode. Finally, from the 
experimental observations of bistability, we know that 
the size of the barrier should be at least 2 kT (Fig. [5]). 

Exploring possible corrections to the calculation, we 
developed a method to include disorder due to the ran- 
domness in the basepair sequence. This disorder intro- 
duces a random intrinsic bend to the DNA, which we are 
able to incorporate by numerically locating the saddle 
point configurations. Intrinsic bends do not significantly 
change the hopping rate, though they do increase the ef- 
fective free energy barrier (Fig. [7]). Both sequence depen- 
dence and uncertainties in the elastic parameters produce 
variations in the rate, but they are not large enough to 
slow the hopping rate by three orders of magnitude to 
the experimental timescale. 

We instead attribute the slowness of the hopping to 
the large bead used to manipulate the DNA, since the 
timescale controlling the bead's motion (oo^ ~ 5 x 10 2 Hz) 
is two orders of magnitude slower than the plectoneme 
nucleation rate (topping ~ 5 x 10 4 Hz). This separa- 
tion of timescales means that the bead moves through 
an effective free energy potential that is set by all pos- 
sible DNA configurations at a given bead position [44] . 
When the bead is near the saddle extension, the energy 
barrier is lower between states with and without a plec- 
toneme (see Fig. [9| , and the microscopic configuration of 
the DNA hops quickly between states with and without 
a plectoneme at a rate faster than topping- But smce 
experiments measure the bead position, this hopping is 
invisible, and we see only the slower hopping of the bead 
(of order 10 Hz), which is set by its own viscous drag. 
The situation is illustrated in Fig. [To] 

If the characteristic rate controlling the bead motion 
were instead made faster than the hopping rate, simi- 
lar experiments could directly measure the plectoneme 
nucleation hopping rate. Some modifications would be 
relatively easy: we estimate that increasing the external 
force F to 5.5 pN (the highest force at which plectonemes 
are observed in the current experiments [20]) would de- 
crease ^hopping by a factor of 3; decreasing the length of 
the DNA to 1 kbp would decrease /popping by a factor 
of about 2; reducing the bead's size to 100 nm would 
increase uj^ by a factor of 2. These modifications would 
bring the two rates closer by about one order of mag- 
nitude. We do not see an obvious way to overcome the 
remaining factor of ten, but leave the challenge open to 
experimentalists. If this challenge can be met, future 
experiments may be able to directly measure the plec- 
toneme nucleation hopping rate, giving useful informa- 
tion about the microscopic dynamics of DNA in water, 
and providing a novel testing ground for transition state 
theory. 
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FIG. 11: (Color online) Local basis vectors. 



Appendix A: The elastic rod model 



FIG. 9: (Color online) Schematic of the double- well in exten- 
sion and plectoneme length coordinates. Integrating over the 
"plectoneme length" dimension would produce the observed 
bistable free energy as a function of extension shown in Fig. [5] 
Note that fixing the bead position at the saddle extension 
produces stable states with and without a plectoneme, with 
a (smaller) free energy barrier between them, producing the 
fast-timescale hopping behavior illustrated in Fig. IlO] 
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FIG. 10: (Color online) Schematic illustrating relevant 
timescales. Since the timescale governing the bead's motion 
through water (2 ms) is much slower than the plectoneme 
nucleation timescale (20 /is), the bead sets the experimen- 
tally measured hopping timescale (100 ms). When the bead 
is located at the saddle point, the free energy barrier to plec- 
toneme formation is lowered (see Fig. [9|, such that plec- 
tonemes nucleate at a rate faster than our calculated /chopping • 
Plectonemes have time to form and disappear many times 
(inset) as the bead moves slowly between the two free energy 
wells. 
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The physical properties of long DNA molecules have 
been found to be well-described by linear elastic the- 
ory (often referred to as the "worm-like chain" model, 
especially in a statistical mechanics context; see, e.g., 
[2TJ[22]). In this formulation, the DNA is modeled as a 
thin elastic rod, and the energy associated with deform- 
ing it from its natural relaxed state is the sum of local 
elastic bending, twisting, and stretching energies. The 
corresponding elastic constants are sensitive to experi- 
mental conditions such as the ionic concentration of the 
surroundings; in our experimental setup, the bend and 
stretch elastic constants B and S can be measured by fit- 
ting force-extension curves, and the (renormalized) twist 
elastic constant C can be measured from the slope of the 
torque as a function of linking number. These values are 
listed in Table [I] [1 . For the low forces in the current ex- 
periment (which are in a biologically-relevant range pQ), 
the stretch elasticity can be safely ignored [45]; we thus 
treat our DNA as an inextensible elastic rod, with energy 
as given in Eq. 0. Furthermore, since all parts of the 
DNA stay sufficiently far from touching each other in the 
saddle state, we also neglect nonlocal repulsive interac- 
tions (which would be required to stabilize plectonemes), 
allowing an analytical description of the saddle state [13] . 



Appendix B: Calculating the energy of a DNA 
configuration 



In our numerical calculations, we approximate the con- 
tinuous elastic rod as a discretized chain of segments, 
each with fixed length d. The orientation of each segment 
is described by the rotations necessary to transform the 
global Cartesian axes onto the local axes of the segment. 
When minimizing the (free) energy, we find it convenient 
to use Euler angles [46] 0, 0, and -0, since any set of Euler 
angles specifies a valid configuration [47]. However, writ- 
ing the energy in terms of differences of Euler angles can 
lead to numerical problems near the singularities at the 
poles. When calculating energies and forces, we therefore 
instead use the full rotation matrix R. R rotates a seg- 
ment lying along the z-axis to its final orientation; i?'s 
columns are thus the two normal unit vectors followed by 



13 



the tangent unit vector (see Fig. 11): 

R M = [ h[ n) 4 n) t» ] (Bl) 

hi = ( cos cos ip — cos sin sin 

— cos sin — cos cos sin -0, 
sin # sin?/;), 

n2 = ( cos cos -0 sin + cos sin -0, (B2) 
cos cos cos -0 — sin sin -0, 

— cos?/; sin 0), 

£ = ( sin sin 0, — cos sin 0, cos 0) . 

There are three degrees of freedom in the "hinge" be- 
tween each segment that determine the local elastic en- 
ergy: the two components of bending fa and fa (along 
h\ and fii, respectively), and the twist T. In terms of the 
rotation matrix 



/\M = (i?( n )) T i^ n+1 



(B3) 



which measures the rotation between adjacent segments, 
mapping the nth segment's axes onto those of the (n + 
l)th, the bends and twist can be written in an explicitly 
rotation-invariant form: to lowest order in the angles (see 
Section D 2 ) , 



/? ■ ni = /3i = (A 23 
/3-n 2 =p 2 = (A 3i 

r = (Ai 2 



A 32 )/2 

A 13 )/2 
A 2 i)/2, 



(B4) 
(B5) 
(B6) 



where the (n) superscripts have been omitted. The above 
forms are useful when the sign of a given component is 
necessary [48 . Otherwise, the following squared forms 
may be used, which have the advantag e of a larger range 
of validity away from zero (see Section |P~2|): 



2 = fi+ft = 2{l 
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Our energy function is then the discretized version of 
Eq. Q: 



E e la 



CT 



(B9) 



or, separating the bend into its two components, 
Mastic = B0 m ■ n { r ] ? + B0 m • n 2 m) ) 2 + CT 
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m ' 



(BIO) 



Appendix C: Transition state calculation details 

1. Changing to the correct coordinates 

So far, we have used Euler angles (0,0,-0) to 
parametrize the DNA configuration. We can easily calcu- 
late energy derivatives with respect to these coordinates. 



However, the hopping rate calculation requires that we do 
our path integral in the same coordinates as the dynam- 
ics, given in Section II C[ which are defined in Cartesian 
space with a local twist [r = (x, y, z, / 0)]. To efficiently 
calculate the correct Hessian, then, we need to convert 
energy derivatives to r space. 

First, we note that there is one less coordinate in Eu- 
ler angle space: this comes from the inextensibility con- 
straint, which r space does not have. Thus we add a 
coordinate A to the Euler angles specifying the length of 
each segment (which will usually be set to a constant d) ; 
we will call this set of coordinates a = (0, 0, -0, A). Also, 
f has TV + 1 elements, each defining the location of one 
end of a segment, while a has N elements, each defining 
the Euler angles and length of each segment. To match 
the number of degrees of freedom, we remove center of 
mass motion and set constant orientation boundary con- 
ditions by fixing the location of the first segment's two 
ends and fixing the last segment's orientation along z: 
r(0) 0, f(l) = dz : r(N + 1) r(N) +dz. This gives a 
total of 4(A^ — 2) degrees of freedom. 

The Jacobian we would like to calculate is 



J rr 



dotrni 

dv n j 



(Cl) 



where m and n label segments and i and j label compo- 
nents. First writing the Euler angles for a given segment 
n in terms of t n = (x, y, z) n+1 - (x, y, z) n , 



= arctan (—t y /t x ) 

= arccos (t z / '^t 2 x + ~t\ + t^j 



(C2) 
(C3) 
(C4) 



We then take derivatives with respect to Cartesian coor- 
dinates to produce the 0, 0, and A rows in the Jacobian. 
A subtlety arises in finding expressions for derivatives of 
the Euler angle ip with respect to Cartesian coordinates. 
We would like the derivative to correspond to rotating the 
adjacent segments to accommodate the change in the lo- 
cation of their connecting ends. Due to the way in which 
Euler angles are defined, this rotation does not in gen- 
eral leave ip unchanged, as one might naively expect. We 
therefore obtain the derivatives of ip by first writing the 
rotation matrix corresponding to infinitesimal motion in 
Cartesian space and then calculating the corresponding 
change in 0. This produces the nonzero terms in the ip 
row of J\ below. 

In the end, we have 



Jmn = ^m,n [Jl(m) + J2] £m,n+l«M m )> ( C5 ) 
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where (including the names of the components for clarity) 
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(C6) 

p = yj A 2 — tj, J2 = 5^,-0, and all the components are 
evaluated at location n. We use this Jacobian to trans- 
form forces with respect to angles a to forces with respect 
to r, which are then used to assemble the Hessian for use 
in calculating the unstable mode and the entropic factors 
for the transition state calculation in Section III Dl 
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2. Other subtleties 

Since derivatives in r space will in general couple to A 
(changing the lengths of segments), we also include an 
extra stretching energy: 



-^stretching ^ ^ ^ (^n • 



(C7) 



FIG. 12: (Color online) Choosing d. Error in the discretized 
energy barrier [comparing to the exact result in Eq. Q] as a 
function of the number of discrete segments N (for L — 740 
nm and other parameters as in Table [I]). Inset: the same data 
on a log-log plot. The green (upper) line has a slope of —2, 
showing convergence proportional to 1/N 2 . We choose N = 
740 (d = 1 nm) as a good trade-off between accuracy and 
required computational resources. 



This avoids problems with extra zero modes correspond- 
ing to changing A. We may choose the stretch elastic 
constant S to match DNA (in which case it should be 
about 1000 pN [H]); but since S is so large compared 
to the other elastic constants, the stretching modes have 
much higher energy and are the same for the straight and 
saddle configurations, canceling in the rate equation [e.g. 
Eq. (21)]. Thus we find that our results are insensitive 
to the exact value of 5, as expected. 

As can be seen by inserting t = (0,0, 1) into Eq. (C6), 



there are singularities in the Jacobian when t points along 
the z-axis. This corresponds to the singularity in the Eu- 
ler angle representation at the poles (when = 0, (j) and 
ip are degenerate) . This is a problem for our formulation 
because our usual boundary conditions hold the ends in 
the z direction. As pointed out in Ref. [23 , a simple 
way to avoid this problem is to rotate the system away 
from the singularity (rotating the direction of the force 
as well). When performing calculations that require the 
Jacobian, we therefore rotate the system about the ?/-axis 
by an angle /3 and modify the external force term in the 
energy from — F cosQ to — F(cos /3cos# + sin/3sin#sin(/>). 
This more general formulation also permits an explicit 
check that our energies and derivatives are rotation in- 
variant. 

The Hessian is constructed by taking numerical deriva- 
tives of forces, which can be calculated analytically. This 
gives the Hessian as a 4(7V — 2) x 4(iV — 2) matrix, which 
is diagonalized to find eigenvalues for the entropic cal- 
culation [Eq. pi])]. At zero disorder, the zero modes 
must first be removed; numerically, we find that the zero 
modes show up conveniently as the two modes with small- 



est eigenvalues, a few orders of magnitude smaller than 
any others. 



Appendix D: Numerical details 

1. Choosing d 

As shown in Fig. [l2j we must be careful to choose our 
discretization length (the length d of each segment) such 
that our energy calculations are sufficiently accurate. We 
can check the accuracy of the discretized energy calcula- 
tion by comparing with the analytical energy barrier in 
Eq. Choosing d = 1 nm (about 3 DNA basepairs) 
produces energy barriers within 0.2/cT of the continuum 
limit (corresponding to 20% changes in the hopping rate) 
with reasonable memory and time expenditure. 



2. Deriving rot at ion- invariant forms for bend and 
twist 

The amount of local bend and twist can be measured 
by differences in the rotation matrices of adjacent seg- 
ments. We would like expressions in terms of the rotation 
matrices that are correct to lowest order in the bend or 
twist angle and that are explicitly rotation invariant. Our 
procedure will be to form rotation invariant terms and, 
writing them in terms of Euler angles and their deriva- 
tives, check what they measure in terms of bend and 
twist. 
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We can first write the bend and twist in terms of 



derivatives of the local basis vectors (see Fig. 11 ) or Euler 
angles [T3] : 



P 2 = 
T 2 = 



= 6 2 sin 2 ( 



r • j2 /. .\2 

(n x n) -t = ( (j> cos + ip J . 



(Dl) 
(D2) 



We then form rotation-invariant terms from rotation 
matrices and compare their Taylor series expansions 
with resp ect to bend and twist angles to Eq. ( |D1[ ) and 
Eq. {mf. First checking Tr[(i^ n+1 ) - RM) t (R^+^ - 
ijWjjTwe find that 



/3 2 



2 V a(3 
3- A n 



c>(n+l) 



*22 



)( 



(n) 

a/3 



E>(n+1) N 



(D3) 



Next we check the dot product of the difference in the 
tangent unit vector t: 



f3 2 = St • St 
= 2(1-, 



(n) 
a3 



(n+l) 



)( 



(n) 
a3 



(n+1)^ 



^33 ) 



(D4) 



Eq. (|D3j) and Eq. (|D4j) produce Eq. <\B7] and Eq. (|B8). 

To find expressions for signed f3 and T, we notice that 
the above use only the diagonal elements of A; we can 
also form rotation invariant terms using the off-diagonal 
elements, producing (with the Levi-Civita symbol e) 

(3r = = (A 23 - A 32 )/2 (D5) 

h = £2 7 sA lS /2 = (A 3 i - A 13 )/2 (D6) 
T = £ 375 A 75 /2 = (A 12 - A 21 )/2. (D7) 

W e ca n see the be nefit of using the squared forms 
[Eq. ( |B7[ ) and Eq. (B8)] by checking their dependence on 
pure bending or twisting rotations. For example, with 
two segments differing only in tw ist, such that ^( n+1 ) = 
ifj^ + a, we find that Eq. (B6) and Eq. (B8) produce, 
respectively [49], T = sin a and T 2 = 2(1 — cos a). Plot- 
ting T 2 for the two cases (Fig. 13) demonstrates that they 



have the same curvature near a — (by design), but us- 
ing the non-squared version in Eq. ( B6 ) leads to a second 



minimum at a = 7r: we find that, especially when includ- 
ing intrinsic bend disorder, this can cause the numerical 
minimizer to allow ip to slip by tt to the next minimum, 
unphysically removing linking number. For this reason, 
we use the squared forms unless otherwise necessary. 



Appendix E: Including intrinsic bend disorder 
1. Setting the disorder size 



As discussed in Section IV, we need to include DNA's 



intrinsic bend disorder to understand the energetics of 
the saddle point barrier crossing. To accomplish this, we 
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FIG. 13: (Color online) Checking bend and twist expressions, 
(top) Two different rotation-invariant approximations to the 
bend or twist angle squared (green, dash-dotted; and red, 
dotted) compared to the actual angle squared (blue, solid). 
Using the non-squared version (green, dash-dotted) leads to a 
smaller range of validity, (bottom) Typical magnitude of bend 
and twist angles for a plectoneme (for F — 2 pN, d — 1 nm). 
Note that the bend and twist angles stay within the region 
where either approximation should be valid. 



shift the zero of the elastic bend energy for each hinge; 
generalizing Eq. ( |B1Q ), 



bend 



= Yd E((^-^)^i m) ) 2 +((^-^)-^ m) ) 2 - 

m 

(El) 

For each z, we choose each of the two components of £ 
from a normal distribution with unit standard deviation. 
We then need to relate to the resulting intrinsic bend 
persistence length P. The persistence length is defined 
by the decay of orientation correlations: 



(t(0) -i(s)} = (cos 0(s)} 



-s/P 



(E2) 



where 0(s) is the angle between segments separate d by 
arclength s [3j. For small s [and thus small Eq. (E2) 



becomes 



1 



(0(s) 2 ) = 1 



s 

P' 



(E3) 



At zero temperature and zero force, the size of (0(s) 2 } is 
set by the intrinsic bends £ m only; we are doing a ran- 
dom walk in two dimensions with one step of root-mean- 
square size \f2o\) taken for every segment of lengt h d . 
Thus (0(s) 2 } =2^a 2 , which when inserted into Eq. ( |E3| 
gives the desired relation between persistence length P 
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FIG. 14: (Color online) Typical sensitivity of saddle energy 
to bending disord er st rength as a function of position [A(s), 
as denned in Eq. (E9)], for L = 740 nm. Peaks indicate po- 
sitions where plectoneme nucleation is energetically favored. 
For DNA, the bending disorder strength D is estimated to 
be of order 0.01 to 0.1 nm -1 / 2 ; thus we expect pinning to 
individual sites (with barriers on the order of kT), but also 
that multiple plectoneme locations will contribute to the fi- 
nal hopping rate in Eq. (32) (with multiple locations having 
energy within about 1 kTj7 



and the size of individual random bends a^: 



(E4) 



We will also define a convenient parameter D control- 
ling disorder strength that is independent of the segment 
length d: 



D . 



(E5) 



such that 



P = D~ 



(E6) 
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FIG. 15: (Color online) The lowest saddle energy as a func- 
tion of disorder strength D. The blue dots are the true saddle 
energies, calculated by numerically zeroing forces on the sad- 
dle configuration. The green line represents the first-order ap- 
proximation given by DA(s*); we see that the approximation 
correctly predicts the scale of the change, but overestimates 
it by many kT at large disorder. 



respect to disorder strength at = 0: 
dE hend 



dE 
da b 



(Tb=0 



dab 
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(E7) 



or, in terms of disorder strength D defined in Eq. (E5) 

B 



dE 
dD 



D=0 



Vd 



(E8) 



We will specifically be interested in the derivative of the 
saddle configuration's energy, which will depend on its 
location s and rotation p. Noting that p will simply ro- 
tate the bend vector /?, we can write down the form of 
the dependence on p: 



2. First order changes in the energy barrier due to 
disorder 



How does disorder change the energy of the saddle 
point? Since the disorder changes only the bending en- 
ergy, we can find the lowest order change fro m the zero 
disorder energy by taking the derivative of Eq. (El ) with 



saddle 



saddle 



dD 



D=0 



-A(s) cos(p-p*(s)) 
(E9) 

for some A(s) and p*(s). A(s) gives the maximum deriva- 
tive (sensitivity to disorder) at position s, and p*(s) is 
the preferred rotation of the saddle that gives the maxi- 
mum (negative) derivative. We can find A(s) and p*(s) 
numerically using the derivative calculated at two values 
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of p separated by w/2: 



M*) = VKaddleM))^ + (E' saddle (s,n/2)Y (E10) 



p*(s) = arctan 



^saddle 



(«,7T/2) 



saddle 



( s ,o) 



(Ell) 



Fig. 14 shows a typical A as a function of s. Fig. 15 com- 
pares the first-order approximation to the saddle energy 
at the location s*, given by DA(s*), to Ca ddie c alculated 
numerically by zeroing forces (see Section E3). We see 
that the approximation correctly predicts the scale of the 
change for the disorder sizes in which we are interested, 
but overestimates the change by many kT for large dis- 
order. We therefore use the first-order approximation 
to find the likely locations of the lowest energy saddle 
points [at the peaks of A(s)], and then zero the forces 
numerically. 

To estimate the typical size of this sensitivity to disor- 
der, we first note that the bend magnitude f or th e saddle 
configuration is [inserting Eqs. ([5| into Eq. (Dl)] 



-'saddle 



2d 

— —seen 



(!)■ 



(E12) 



Approximating the function sech (x) as 1 in the range 



-2 < x < 2 and elsewhere, Eq. (E8) becomes 



^saddle 



dD 



D=0 



B 2d 



U/d 
n=0 



(E13) 



for randomly distributed £ with unit variance, which pro- 
duces 



/ / ^saddle \2v 4£? 

U dD } yft 



(E14) 



Inserting B = (43 nm)/cT and £ ~ 10 nm gives a typical 
sensitivity of about (50 nm 1//2 )/cT, agreeing with the scale 



found in the full calculation, as shown in Fig. 14 



3. Finding saddle points 

With large disorder, the saddle points must be found 
numerically. We start by estimating the set of saddle lo- 
cations {s*} and rotat ions {p*} using first order pertur- 
bation theory (Section |E 2| ). We find the local maxima of 
A(s) using a one-dimensional local search method start- 
ing from a set of points spaced by the saddle configuration 
length scale £. This gives {s*}, from which {p*} can be 
found using Eq. (Ell ). For each s* and corresponding p*, 
we create a zero-disorder saddle configuration [Eqs. ([5|], 
and then use this as the starting point for a multidimen- 
sional equation solver (scipy . f solve) that numerically 
locates the saddle with disorder by finding solutions with 
zero net force on each segment. 
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the Oseen-Burgers tensor [25]) that is modified at short 
distances such that it becomes positive definite, produc- 
ing stable dynamics |26| . 

Choosing these units for our variables makes the 
path integral partition function unitless: Z 

C T r dx n dy n dz n dt/j n —E({x,y,z,ip})/kT 

J 1 lr7 '» 6 



I 3 



The rate will be slower if the system must traverse ; 
narrow "pass" at the saddle point. 
Using M diagonal [Eq. (13)] prod uces a comparable unsta- 



ble mode rate of Xp/^Ti — 5.6 X 10 Hz. We use the more 
physical M R 
calculations. 
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